module quickbeam
  USE COSP_KINDS,           ONLY: wp
  USE MOD_COSP_CONFIG,      ONLY: DBZE_BINS,DBZE_MIN,DBZE_MAX,CFAD_ZE_MIN,CFAD_ZE_WIDTH, &
                                  R_UNDEF,cloudsat_histRef,use_vgrid,vgrid_zl,vgrid_zu
  USE MOD_COSP_STATS,       ONLY: COSP_LIDAR_ONLY_CLOUD,hist1D,COSP_CHANGE_VERTICAL_GRID
  implicit none

  integer,parameter :: &
       maxhclass     = 20,  & ! Qucikbeam maximum number of hydrometeor classes.
       nRe_types     = 550, & ! Quickbeam maximum number or Re size bins allowed in N and Z_scaled look up table.
       nd            = 85,  & ! Qucikbeam number of discrete particles used in construction DSDs.
       mt_ntt        = 39,  & ! Quickbeam number of temperatures in mie LUT.
       Re_BIN_LENGTH = 10,  & ! Quickbeam minimum Re interval in scale LUTs  
       Re_MAX_BIN    = 250    ! Quickbeam maximum Re interval in scale LUTs
  real(wp),parameter :: &
       dmin          = 0.1, & ! Quickbeam minimum size of discrete particle
       dmax          = 10000. ! Quickbeam maximum size of discrete particle
  
  logical :: radar_at_layer_one   ! If true radar is assume to be at the edge 
                                  ! of the first layer, if the first layer is the
                                  ! surface than a ground-based radar.   If the
                                  ! first layer is the top-of-atmosphere, then
                                  ! a space borne radar.

  ! ##############################################################################################
  type radar_cfg
     ! Radar properties
     real(wp) :: freq,k2
     integer  :: nhclass                           ! number of hydrometeor classes in use
     integer  :: use_gas_abs, do_ray
 
     ! Variables used to store Z scale factors
     character(len=240)                             :: scale_LUT_file_name
     logical                                        :: load_scale_LUTs, update_scale_LUTs
     logical, dimension(maxhclass,nRe_types)        :: N_scale_flag
     logical, dimension(maxhclass,mt_ntt,nRe_types) :: Z_scale_flag,Z_scale_added_flag
     real(wp),dimension(maxhclass,mt_ntt,nRe_types) :: Ze_scaled,Zr_scaled,kr_scaled
     real(wp),dimension(maxhclass,nd,nRe_types)     :: fc, rho_eff
     real(wp),dimension(Re_MAX_BIN)                 :: base_list,step_list

  end type radar_cfg

contains
  ! ##############################################################################################
  !                                     QUICKBEAM_PIXEL
  ! ##############################################################################################
  subroutine quickbeam_pixel(rcfg,nprof,ngate,undef,hgt_matrix,z_vol,kr_vol,g_vol, &
                             a_to_vol,g_to_vol,dBZe,Ze_non,Ze_ray)

    ! INPUTS
    type(radar_cfg),intent(inout) :: &
         rcfg             ! Derived type for radar simulator setup
    integer,intent(in) :: &
         nprof,         & ! Number of hydrometeor profiles
         ngate            ! Number of vertical layers
    real(wp),intent(in) :: &
         undef            ! Missing data value
    real(wp),intent(in),dimension(nprof,ngate) :: &
         hgt_matrix,    & ! Height of hydrometeors (km)
         z_vol,         & ! Effective reflectivity factor (mm^6/m^3)
         kr_vol,        & ! Attenuation coefficient hydro (dB/km)
         g_vol            ! Attenuation coefficient gases (dB/km)

    ! OUTPUTS
    real(wp), intent(out),dimension(nprof,ngate) :: &
         Ze_non,        & ! Radar reflectivity without attenuation (dBZ)
         Ze_ray,        & ! Rayleigh reflectivity (dBZ)
         g_to_vol,      & ! Gaseous atteunation, radar to vol (dB)
         a_to_vol,      & ! Hydromets attenuation, radar to vol (dB)
         dBZe             ! Effective radar reflectivity factor (dBZ)

    ! LOCAL VARIABLES
    integer :: k,pr,start_gate,end_gate,d_gate
    real(wp),dimension(nprof,ngate) :: &
         z_ray      ! Reflectivity factor, Rayleigh only (mm^6/m^3)

    ! Load scaling matricies from disk -- but only the first time this subroutine is called
    if(rcfg%load_scale_LUTs) then
       call load_scale_LUTs(rcfg)
       rcfg%load_scale_LUTs=.false.
       rcfg%Z_scale_added_flag = .false. ! will be set true if scaling Look Up Tables are modified during run
    endif

    ! Initialization
    g_to_vol = 0._wp
    a_to_vol = 0._wp

    ! Loop over each range gate (ngate) ... starting with layer closest to the radar !
    if(radar_at_layer_one) then
       start_gate = 1
       end_gate   = ngate
       d_gate     = 1
    else
       start_gate = ngate
       end_gate   = 1
       d_gate     = -1
    endif
    do k=start_gate,end_gate,d_gate
       ! Loop over each profile (nprof)
       do pr=1,nprof
          ! Attenuation due to hydrometeors between radar and volume
          
          ! NOTE old scheme integrates attenuation only for the layers ABOVE
          ! the current layer ... i.e. 1 to k-1 rather than 1 to k ...
          ! which may be a problem.   ROJ
          ! in the new scheme I assign half the attenuation to the current layer
          if(d_gate==1) then
             ! dheight calcuations assumes hgt_matrix points are the cell mid-points.
             if (k>2) then
                ! add to previous value to half of above layer + half of current layer
                a_to_vol(pr,k)=  a_to_vol(pr,k-1) + &
                     (kr_vol(pr,k-1)+kr_vol(pr,k))*(hgt_matrix(pr,k-1)-hgt_matrix(pr,k))
             else
                a_to_vol(pr,k)=  kr_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k+1))
             endif
          else   ! d_gate==-1
             if(k<ngate) then
                ! Add to previous value half of above layer + half of current layer
                a_to_vol(pr,k) = a_to_vol(pr,k+1) + &
                     (kr_vol(pr,k+1)+kr_vol(pr,k))*(hgt_matrix(pr,k+1)-hgt_matrix(pr,k))
             else
                a_to_vol(pr,k)= kr_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k-1))
             endif
          endif
       enddo   ! End loop over pr (profile)
    enddo ! End loop of k (range gate)
          
    do k=start_gate,end_gate,d_gate
       ! Loop over each profile (nprof)
       do pr=1,nprof
          ! Attenuation due to gaseous absorption between radar and volume
          if ( (rcfg%use_gas_abs == 1) .or. ((rcfg%use_gas_abs == 2) .and. (pr == 1)) ) then
             !DS MOVED UP  g_vol(pr,k) = gases(p_matrix(pr,k),t_kelvin,rh_matrix(pr,k),rcfg%freq)
             if (d_gate==1) then
                if (k>1) then
                   ! Add to previous value to half of above layer + half of current layer
                   g_to_vol(pr,k) =  g_to_vol(pr,k-1) + &
                        0.5*(g_vol(pr,k-1)+g_vol(pr,k))*(hgt_matrix(pr,k-1)-hgt_matrix(pr,k))
                else
                   g_to_vol(pr,k)=  0.5_wp*g_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k+1))
                endif
             else   ! d_gate==-1
                if (k<ngate) then
                   ! Add to previous value to half of above layer + half of current layer
                   g_to_vol(pr,k) = g_to_vol(pr,k+1) + &
                        0.5_wp*(g_vol(pr,k+1)+g_vol(pr,k))*(hgt_matrix(pr,k+1)-hgt_matrix(pr,k))
                else
                   g_to_vol(pr,k)= 0.5_wp*g_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k-1))
                endif
             endif
          elseif(rcfg%use_gas_abs == 2) then
             ! Using value calculated for the first column
             g_to_vol(pr,k) = g_to_vol(1,k)
          elseif (rcfg%use_gas_abs == 0) then
             g_to_vol(pr,k) = 0._wp
          endif
       enddo   ! End loop over pr (profile)
    enddo ! End loop of k (range gate)
    
    ! Compute Rayleigh reflectivity, and full, attenuated reflectivity
    if(rcfg%do_ray == 1) then
       Ze_ray(1:nprof,1:ngate) = merge(10._wp*log10(z_ray(1:nprof,1:ngate)), 1._wp*undef, z_ray(1:nprof,1:ngate) > 0._wp)
    else 
      Ze_ray(1:nprof,1:ngate) = undef
    end if

    where(z_vol(1:nprof,1:ngate) > 0._wp) 
      Ze_non(1:nprof,1:ngate) = 10._wp*log10(z_vol(1:nprof,1:ngate))
      dBZe(1:nprof,1:ngate) = Ze_non(1:nprof,1:ngate)-a_to_vol(1:nprof,1:ngate)-g_to_vol(1:nprof,1:ngate)
    elsewhere
      dBZe(1:nprof,1:ngate) = undef
      Ze_non(1:nprof,1:ngate) = undef
    end where 

    ! Save any updates made 
    if (rcfg%update_scale_LUTs) call save_scale_LUTs(rcfg)
 
  end subroutine quickbeam_pixel
  ! ##############################################################################################
  !                                     QUICKBEAM_GRIDMEAN
  ! ##############################################################################################
  subroutine quickbeam_gridbox(npoints, ncolumns, nlevels,llm, tmp, pnorm, pnorm_perp,   &
                               pmol, Ze_tot,zlev, zlev_half, cfad_ze,                    &
                               lidar_only_freq_cloud,radar_lidar_tcc)
    ! Inputs
    integer,intent(in) :: &
         npoints,    & ! Number of horizontal grid points
         ncolumns,   & ! Number of subcolumns
         nlevels,    & ! Number of vertical layers in OLD grid
         llm           ! NUmber of vertical layers in NEW grid
    real(wp),intent(in),dimension(npoints,ncolumns,Nlevels) :: &
         pnorm,      & ! Lidar ATB
         pnorm_perp, & ! Lidar perpendicular ATB
         Ze_tot        ! 
    real(wp),intent(in),dimension(npoints,Nlevels) :: &
         pmol,       & ! Molecular ATB
         tmp,        & ! Temperature at each levels
         zlev,       & ! Model full levels
         zlev_half     ! Model half levels
         
    ! Outputs
    real(wp),intent(inout),dimension(npoints,DBZE_BINS,llm) :: &
         cfad_ze    !
    real(wp),intent(inout),dimension(npoints,llm) :: &
         lidar_only_freq_cloud
    real(wp),intent(inout),dimension(npoints) :: &
         radar_lidar_tcc

    ! Local variables
    integer :: i,j
    real(wp),dimension(npoints,ncolumns,llm) :: ze_totFlip
    
    if (use_vgrid) then
       ! Regrid in the vertical
       call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,zlev,zlev_half,Ze_tot, &
            llm,vgrid_zl,vgrid_zu,Ze_totFlip,log_units=.true.)

       ! Effective reflectivity histogram
       do i=1,Npoints
          do j=1,llm
             cfad_ze(i,:,j) = hist1D(Ncolumns,Ze_totFlip(i,:,j),DBZE_BINS,cloudsat_histRef)
          enddo
       enddo
       where(cfad_ze .eq. 0)       cfad_ze = R_UNDEF
       where(cfad_ze .ne. R_UNDEF) cfad_ze = cfad_ze/Ncolumns

    else
       ! Effective reflectivity histogram
       do i=1,Npoints
          do j=1,llm
             cfad_ze(i,:,j) = hist1D(Ncolumns,Ze_tot(i,:,j),DBZE_BINS,cloudsat_histRef)
          enddo
       enddo
       where(cfad_ze .eq. 0)       cfad_ze = R_UNDEF
       where(cfad_ze .ne. R_UNDEF) cfad_ze = cfad_ze/Ncolumns
    endif   

  end subroutine quickbeam_gridbox

  
  ! ##############################################################################################
  ! ##############################################################################################

  ! ##############################################################################################

  ! ##############################################################################################
  subroutine load_scale_LUTs(rcfg)
    
    type(radar_cfg), intent(inout) :: rcfg
    logical                        :: LUT_file_exists
    integer                        :: i,j,k,ind
    
    ! Load scale LUT from file 
    inquire(file=trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat', &
         exist=LUT_file_exists)
    
    if(.not.LUT_file_exists) then  
       write(*,*) '*************************************************'
       write(*,*) 'Warning: Could NOT FIND radar LUT file: ', &
            trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'        
       write(*,*) 'Will calculated LUT values as needed'
       write(*,*) '*************************************************'
       return
    else
       OPEN(unit=12,file=trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat',&
            form='unformatted', &
            err= 89, &
            access='DIRECT',&
            recl=28)
       write(*,*) 'Loading radar LUT file: ', &
            trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
       
       do i=1,maxhclass
          do j=1,mt_ntt
             do k=1,nRe_types
                ind = i+(j-1)*maxhclass+(k-1)*(nRe_types*mt_ntt)
                read(12,rec=ind) rcfg%Z_scale_flag(i,j,k), &
                     rcfg%Ze_scaled(i,j,k), &
                     rcfg%Zr_scaled(i,j,k), &
                     rcfg%kr_scaled(i,j,k)
             enddo
          enddo
       enddo
       close(unit=12)
       return 
    endif
    
89  write(*,*) 'Error: Found but could NOT READ radar LUT file: ', &
         trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
    stop
  end subroutine load_scale_LUTs
  
  ! ##############################################################################################
  ! ##############################################################################################
  subroutine save_scale_LUTs(rcfg)
    type(radar_cfg), intent(inout) :: rcfg
    logical                        :: LUT_file_exists
    integer                        :: i,j,k,ind
    
    inquire(file=trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat', &
         exist=LUT_file_exists)
    
    OPEN(unit=12,file=trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat',&
         form='unformatted',err= 99,access='DIRECT',recl=28)
    
    write(*,*) 'Creating or Updating radar LUT file: ', &
         trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
    
    do i=1,maxhclass
       do j=1,mt_ntt
          do k=1,nRe_types
             ind = i+(j-1)*maxhclass+(k-1)*(nRe_types*mt_ntt)
             if(.not.LUT_file_exists .or. rcfg%Z_scale_added_flag(i,j,k)) then
                rcfg%Z_scale_added_flag(i,j,k)=.false.
                write(12,rec=ind) rcfg%Z_scale_flag(i,j,k), &
                     rcfg%Ze_scaled(i,j,k), &
                     rcfg%Zr_scaled(i,j,k), &
                     rcfg%kr_scaled(i,j,k)
             endif
          enddo
       enddo
    enddo
    close(unit=12)
    return 
    
99  write(*,*) 'Error: Unable to create/update radar LUT file: ', &
         trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
    return  
    
  end subroutine save_scale_LUTs
  ! ##############################################################################################
  ! ##############################################################################################
  subroutine quickbeam_init()

    
  end subroutine quickBeam_init
  ! ##############################################################################################
  ! ##############################################################################################


end module quickbeam


